% Analysis Reduction Error 
clc;
clear all;
close all;
%% function analysis_red_err_alpha(lower,upper,m,n,r,k, alpha)
m = 16; n = 16; r = 16; k = 1;
lower = 0;
upper = 1;
ulp = 2^(-r);
step_a = 2^(-m);
step_b = 2^(-n);

num = (round((upper-lower)/step_a))*(round((upper-lower)/step_b));

MAE = zeros(1, 6);
MXE = zeros(1, 6);

%%
id = 0;
for a = lower:step_a:upper
    for b = lower:step_b:upper        
        %answers true, uniform constant and opt. constant
        id = id + 1;
        fprintf('[%d/%d] = %2.2f\n', id, num, double(100 * id/num));
        ans_true = a * b;
        
        ans_rnd = rnd_near(ans_true, r);
        err_rnd = abs(ans_rnd - ans_true);
        MAE(1) = MAE(1) + err_rnd;
        if (err_rnd > MXE(1))
            MXE(1) = err_rnd;
        end
        
        ans_trunc = trunc_mult(a, b,0.0,m, n, r, k);
        err_trunc = abs(ans_trunc - ans_true);
        MAE(2) = MAE(2) + err_trunc;
        if (err_trunc > MXE(2))
            MXE(2) = err_trunc;
        end
        [const_schulte, const_lim] = ccm_const(m, n, r, k);
        ans_lim = trunc_mult(a, b,const_lim, m, n, r, k);
        err_lim = abs(ans_lim - ans_true);
        MAE(3) = MAE(3) + err_lim;
        if (err_lim > MXE(3))
            MXE(3) = err_lim;
        end        
        ans_schulte = trunc_mult(a, b, const_schulte, m, n, r, k);
        err_schulte = abs(ans_schulte - ans_true);
        MAE(4) = MAE(4) + err_schulte;
        if (err_schulte > MXE(4))
            MXE(4) = err_schulte;
        end          
        
        ans_king = vcm(a, b, m, n, r, k);
        err_king = abs(ans_king - ans_true);
        MAE(5) = MAE(5) + err_king;
        if (err_king > MXE(5))
            MXE(5) = err_king;
        end         
        ans_stine = hcm(a, b, m, n, r, k, 0.8);
        err_stine = abs(ans_stine - ans_true);
        MAE(6) = MAE(6) + err_stine;
        if (err_stine > MXE(6))
            MXE(6) = err_stine;
        end        
    end
end
MAE = MAE./num;
MAE = MAE./ulp;
MXE = MXE./ulp;

% %% Display results
% figure(1);
% plot(alpha,abs(err_uniform_alpha), alpha, abs(err_opt_alpha)); 
% legend('Uniform', 'Optimal');
